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Integral  equation  calculations  were  performed  on  two  model  supercritical  mixtures:  naphthalene  in  carbon 
dioxide  (non-volatile  solute;  attractive  mixture)  and  neon  in  xenon  (volatile  solute;  repulsive  mixture).  Both 
systems  were  studied  at  high  dilution,  at  supercritical  temperature,  and  over  a  broad  range  of  densities.  The 
attractive  system  exhibited  significant  short-ranged  solvent  enrichment  around  the  solute.  The  difference 
between  the  solvent  density  averaged  over  three  solvation  shells,  and  the  bulk  solvent  density  was  found  to 
be  more  pronounced  between  50  and  80%  of  the  solvent's  critical  density.  In  contrast,  the  repulsive  system 
exhibited  local  solvent  depletion  near  the  solvent's  critical  point  The  extent  of  the  solvation  region 
responsible  for  determining  the  solute's  chemical  potential  varies  between  three  and  live  solvent  diameters 
for  the  attractive  mixture,  with  the  maximum  occurring  near  the  solvent's  critical  density.  At  60%  of  the 
solvent's  critical  density,  the  local  concentration  of  naphthalene  molecules  around  each  other  is  almost  ten 
times  higher  than  in  the  bulk.  Dynamic  simulations  of  pyrene  in  carbon  dioxide  revealed  insignificant 
solute-solute  aggregation.  The  mean  lifetime  of  transient  solute  dimers,  roughly  0.8  picoseconds,  was  found 
not  to  change  even  when  the  bulk  density  was  varied  by  a  factor  of  four. 
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1.  Introduction 


A  fluid  is  said  to  be  supercritical  when  its  temperature  and  pressure  are 
simultaneously  higher  than  their  critical  point  values  (McHugh  and  Krukonis,  1986).  In 
practice,  though,  the  term  supercritical  is  used  to  describe  fluids  within  the  relatively 
narrow  range  of  temperatures  and  pressures  (1  <  T/Tc  <1.1;  1  <  P/Pc  <  2),  where 
subscript  c  denotes  critical  point  values  (Brennecke  and  Eckert,  1989).  It  is  within  this 
region  that  most  property  changes  involved  in  the  transition  from  a  dilute  gas  to  a  dense 
fluid  occur.  Consequently,  the  thermophysical  properties  of  supercritical  fluids  exhibit  high 
rates  of  change  with  respect  to  temperature  and  pressure.  Not  all  properties  change  at  the 
same  rate;  hence  the  combination  of  thermophysical  properties  of  supercritical  fluids  is 
quite  unique.  For  example,  supercritical  fluids  have  liquid-like  densities,  kinematic 
viscosities  that  are  lower  than  those  of  liquid  metals,  compressibilities  that  can  be  arbitrarily 
higher  than  for  an  ideal  gas,  and  viscosities  that  are  intermediate  between  gas-  and  liquid¬ 
like  (Debenedetti  and  Reid,  1986). 


Dilute  mixtures  in  which  the  solvent  is  supercritical  can  be  classified  into  three 
categories:  attractive,  weakly  attractive,  and  repulsive  (Petsche  and  Debenedetti,  1991).  The 
first  are  of  technological  interest,  as  nearly  all  applications  of  supercritical  fluids  involve 
attractive  solutions.  In  such  a  class  of  binary  mixture,  there  is  long-ranged  solvent 
enrichment  around  the  solute  molecules  in  the  vicinity  of  the  solvent’s  critical  point, 
solubility  enhancement,  and  the  solute's  partial  molar  volume  and  enthalpy  are  large  and 
negative,  diverging  to  -  <*>  at  the  solvent's  critical  point  in  the  limit  of  infinite  dilution. 
Attractive  behavior  occurs  generally  when  the  solute  is  larger  than  the  solvent  and  has  a 


correspondingly  larger  characteristic  interaction  energy.  In  practical  terms,  this  generally 
means  a  low-volatility  solute  dissolved  in  a  supercritical  solvent.  Repulsive  solutions  are 
mainly  of  theoretical  interest,  although  experiments  with  repulsive  supercritical  solutions 
have  been  reported  (Biggerstaff  and  Wood,  1988  a,b).  In  a  repulsive  mixture,  there  is 
long-ranged  solvent  depletion  around  the  solute  molecules  in  the  neighborhood  of  the 


solvent's  critical  point,  and  the  solute's  partial  molar  volume  and  enthalpy  are  large  and  'pJJJT 
positive,  diverging  to  +  °°  at  the  solvent's  critical  point  in  the  limit  of  infinite  dilution,  i 
Repulsive  behavior  occurs  when  the  solvent  is  larger  than  the  solute  and  has  a 
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correspondingly  larger  characteristic  interaction  energy.  The  microstructure  around  solute  u*_ 
molecules  tends  to  be  solvent-rich  in  the  attractive  case,  and  solvent-lean  for  repulsive 
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A  large  body  of  experimental  (Kim  and  Johnston,  1987a, b;  Johnston  et  al„  1989; 
Brennecke  and  Eckert,  1989;  Brennecke  et  al.,  1990  a,b;  Betts  et  al.,  1992  a,b;  Sun  et  al., 

1992) ,  theoretical  (Wu  et  al.,  1990;  Munoz  and  Chimowitz,  1992;  Tom  and  Debenedetti, 

1993) ,  and  computadonal  (Petsche  and  Debenedetti,  1989;  Knutson  et  al.,  1992;  O’Brien  et 
al.,  1993)  evidence  suggests  that  in  dilute  supercritical  mixtures  of  practical  interest  the 
local  environment  surrounding  solute  molecules  (microstructure)  differs  appreciably  from 
the  bulk.  Thus,  an  improved  understanding  of  the  structure  and  dynamics  of  the  local 
environment  around  solute  species  could  well  lead  to  the  tailoring  of  specific  solvent 
environments  for  targeted  separations  or  reactions. 


2.  Objective 

Potential  applicadons  of  supercritical  fluids  range  from  their  use  as  highly  selective 
solvents  for  difficult  separations  (Brennecke  and  Eckert,  1989),  to  materials  processing  and 
particle  formation  (Tom  and  Debenedetti,  1991).  The  actual  implementation  of  some  of 
these  promising  technologies  has  been  hampered  by  a  lack  of  understanding  of  key  aspects 
of  the  behavior  of  fluids  and  their  mixtures  at  supercritical  conditions,  such  as  molecular 
interactions  and  solvation  mechanisms,  transport  properties  and  reaction  kinetics,  and 
mechanisms  of  nucleation  and  growth  leading  to  particle  formation.  The  objective  of  this 
research  is  to  gain  an  understanding  of  the  microscopic  structure  and  dynamical  behavior  of 
supercritical  mixtures,  using  integral  equation  calculations  and  molecular  dynamics 
computer  simulations.  The  quantities  of  interest  include  molecular  distribution  functions, 
local  density  augmentation  and  depletions,  solute  rotational  dynamics,  solvation 
mechanisms  and  dynamics,  dynamics  of  solute-solute  collisions,  solute  translational 
diffusion.  In  this  report,  we  describe  integral  equation  calculations  and  molecular  dynamics 
simulations  of  model  attractive  and  repulsive  supercritical  Lennard-Jones  mixtures. 


3.  Status  of  the  Research  Effort 
Integral  Equation  Calculations 

Integral  equations  are  useful  for  studying  dilute  supercritical  mixtures  because 
arbitrarily  low  solute  mole  fractions  can  be  studied  at  no  additional  computational  cost.  The 
technique  only  yields  structural,  thermodynamic  information.  Two  highly  dilute  binary 
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supercritical  systems  were  studied:  one  is  attractive  and  one  repulsive.  The  spherically 
symmetric  Lennard-Jones  potential  was  used  in  both  cases.  Potential  parameters  are  listed 
in  Tables  I  and  II.  These  mixtures  have  been  called  naphthalene  in  carbon  dioxide 
[attractive  case;  see  Cochran  and  Lee  (1989)]  and  neon  in  xenon  [repulsive  case;  see 
Petsche  and  Debenedetti  (1989)].  We  use  this  terminology,  although  it  should  be 
understood  that  except  for  yielding  reasonable  estimates  for  the  pure-component  critical 
temperatures  and  densities  with  the  parameters  listed  in  Tables  I  and  II,  the  Lennard-Jones 
potential  is  a  rudimentary  representation  of  the  actual  interaction  potentials,  especially  in  the 
attractive  case.  The  basic  features  of  the  attractive  and  repulsive  behavior  to  be  discussed 
below,  however,  are  quite  general,  and  hence  independent  of  the  details  of  the 
intermolecular  potentials. 

Table  I.  Lennard-Jones  Parameters  for  Naphthalene  in  CO2 


Naphthalene 
in  CO2 

9 

Oij  (A) 

Oij/022 

eij/k  (K) 

£ij/E22 

solute-solute 

11 

6.199 

1.634 

554.4 

2.458 

solvent-solvent 

22 

3.794 

1.000 

225.5 

1.000 

solute-solvent 

12 

4.996 

1.317 

353.4 

1.567 

Table  II.  Lennard-Jones  Parameters  for  Neon  in  Xenon 


Neon 

in  Xenon 

U 

Cjj  (A) 

Oij/^22 

Eij/k  (K) 

Eij/622 

solute-solute 

11 

2.82 

0.697 

32.8 

0.142 

solvent-solvent 

22 

4.047 

1.000 

231.0 

1.000 

solute-solvent 

12 

3.433 

0.848 

87.0 

0.377 
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Dimensionless  temperatures,  densities,  and  distances  are  denoted  by  a  superscript  *,  and 

are  given  by  T*  =  kT/£2,  p*  =  and  r*  =  r/C2,  respectively  (1  =  solute;  2  =  solvent).  A 

solute  mole  fraction  of  10'9  was  used  in  all  the  calculations  to  approximate  infinitely  dilute 

conditions.  Attractive  mixture  calculations  were  done  at  T*  =  1.4,  and  0.1  <  p*  <  0.5. 

Repulsive  mixture  calculations  were  done  at  T*  =  1.4  and  0.1  <  p*  <  0.7.  Using  the  best 

*  _* 

estimate  of  the  critical  point  for  a  Lennard-Jones  fluid  ( pc  -  0.31;  Tc  =  1.31;  Smit  et  al., 
1989)  this  corresponds  to  a  reduced  temperature,  T/Tc  =  1 .07  (attractive  and  repulsive 
cases),  and  reduced  density  ranges  (p/pc)  of  0.32  to  1.61  (attractive  case)  and  0.32  to  2.26 
(repulsive  case). 

Figures  1  and  2  show  the  solute-solute  (11)  and  solute-solvent  (12)  pair  correlation 
functions  for  the  attractive  mixture.  The  first  peak  of  the  1 1  correlation  occurs  at  r*  =  1.8, 
and  represents  the  direct  contact  between  naphthalene  molecules.  This  distance  corresponds 
to  the  well  of  the  1 1  potential,  which  occurs  at  a  separation  of  1.83.  The  maximum  peak 
occurs  at  a  reduced  density  (p/pc)  of  0.58.  A  second  feature  of  the  solute-solute  correlation 
function  is  a  broad  shoulder  occurring  approximately  between  r*  =  2.3  and  3  at  sub-critical 
densities.  For  densities  above  p*  =  0.25,  the  shoulder  becomes  a  second  peak,  whose 
location  (r*  =  2.76)  suggests  a  contribution  from  naphthalene-C02-naphthalene  triplets. 
The  minimum  between  the  first  two  peaks  decreases  with  density:  it  falls  below  unity  at  p* 
=  0.4.  These  calculations  provide  further  evidence  of  high  solute-solute  peaks  in  attractive 
mixtures  (Chialvo  and  Debenedetti,  1992),  especially  at  sub-critical  densities.  The  solute- 
solvent  (12)  correlation  (Figure  2)  shows  a  first  peak  at  r*  =  1.44,  representing  the  direct 
contact  between  naphthalene  and  CO2  molecules.  Note  that,  contrary  to  the  1 1  function,  the 
first  peak  attains  its  highest  value  (gj2  =  3.07)  at  the  lowest  density  studied  here  (p*  =  0.1, 
vs.  p*  =  0.18  for  the  gj  1  maximum  peak). 

A  convenient  way  of  quantifying  the  local-bulk  asymmetry  of  the  distribution  of 
solvent  molecules  is  by  integrating  the  solute-solvent  pair  correlation  function  and  defining 
a  radially-dependent  local  density  within  a  sphere  of  radius  r  around  a  central  solute 
molecule. 


<  N2 1  (r )  > =p2  J 4  nf '2  gn  (r'  )dr‘ 
0 


(1) 
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Piocal  (r )' 


<Nn(r)> 
— - - - - — 

i*  r3_f  OjY 

3  l  2  J 


(2) 


where  <N  i2(r)>  is  the  average  number  of  solvent  molecules  within  a  sphere  of  radius  r 
centered  around  a  solute  molecule,  and  piocal  is  the  average  solvent  density  within  the  same 
sphere,  but  excluding  the  solute's  effective  volume.  Figure  3  shows  the  radial  and  density 
dependence  of  piocal  for  the  attractive  mixture,  normalized  by  the  bulk  density  (so  that  all 
curves  decay  to  unity  at  sufficiently  large  distances).  The  radius  of  the  partially  hollow 
sphere  within  which  density  augmentation  is  largest  decreases  monotonically  with  density, 
from  r*  =  1.84  at  p*  =0.1  to  r*  =  1.72  at  p*  =  0.5.  The  maximum  density  enhancement 
always  occurs  within  the  first  solvation  shell.  However,  significant  density  augmentation 
persists  even  up  to  three  solvation  shells.  At  low  and  near-critical  densities,  the 
enhancement  persists  over  large  distances  away  from  the  solute  molecule  (for  example,  at 
p*  =  0.2,  piocal/Pbulk  =  1  045  at  r*  =  8).  However,  ai  higher  densities,  the  enhancement 
decays  much  more  rapidly  (e.g.,  at  p*  =  0.5,  piocal/Pbulk  =  101  at  r*  =  5).  This  is  due  to  a 
cancellation  between  peaks  and  valleys  in  the  correlation  function  as  liquid-like  densities  are 
approached.  Thus,  we  see  an  augmentation  in  the  solvent's  local  density  that  is  particularly 
pronounced  at  subcritical  densities  and  decreases  rapidly  at  liquid-like  densities. 

Corresponding  calculations  are  shown  in  Figures  4,  5  for  the  repulsive  mixture. 
The  solute-solvent  pair  correlation  function  shows  a  marked  difference  with  respect  to  the 
corresponding  attractive  quantity  (Figure  2).  In  the  attractive  case  we  see  a  pronounced  first 
peak  that  decreases  with  density,  and  a  shallow  first  valley.  In  the  repulsive  case,  we  see  a 
mild  first  peak  which  increases  with  density,  and  a  deep  first  valley.  The  density  and  radial 
dependence  of  the  local  density  [as  defined  in  (2)]  is  shown  in  Figure  5.  Note  that  at  p*  = 
0.1,  piocal/Pbulk  approaches  i.O  monotonically  from  below.  Above  this  density,  the 
monotonic  rise  evolves  gradually,  first  into  a  curve  with  an  inflection  point  at  r*  =  1.75, 
and  then  into  a  curve  with  maxima  and  minima.  Minima  and  inflection  points  occur 
between  r*  =  1.6  and  1.8  (within  the  second  solvation  shell),  with  the  deepest  minimum 
occurring  at  a  density  p*  =  0.3.  Thus,  unlike  the  case  of  the  attractive  mixtures  (where  the 
most  pronounced  local  density  enhancements  are  found  at  substantially  subcritical 
densities),  repulsive  mixtures  show  the  greatest  density  depletion  close  to  the  critical 
density.  Correlation  holes  such  as  the  one  shown  in  Figure  5  are  unusual.  They  occur  here 
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not  because  of  packing  constraints  (the  usual  cause  of  correlation  holes  in  branched  and 
chain  molecules),  but  because  of  competitive  attractions  (£22>>el2)  at  moderate  density. 
Note  the  rapid  disappearance  of  repulsive  behavior  at  high  density. 


Integral  equations  provide  a  link  between  microstructure  and  thermodynamics 
because  they  allow  the  calculation  of  chemical  potentials.  Specifically,  we  are  interested  in 
the  solute's  fugacity  coefficient,  <j>,,  defined  by 


H\=kTln 


y\<hpA] 

kT 


(3) 


where  pi  is  the  solute's  chemical  potential,  yj,  its  mole  fraction  in  the  fluid  phase,  P  is  the 
pressure,  Ai  is  the  solute's  deBroglie  wavelength,  and  k  is  Boltzmann’s  constant.  The 
fugacity  coefficient,  in  other  words,  is  the  proportionality  constant  between  partial  pressure 
and  fugacity 


fl  =  y\P<P\ 


(4) 


This  quantity  is  of  fundamental  importance  for  solubility  calculations.  In  particular,  it  is 
related  to  yi,  the  equilibrium  mole  fraction  of  a  pure,  incompressible,  non-volatile  solute  in 
a  fluid  phase  by 


0-1  _  y\{p  1  Pvap) 

exP[pv  solid  1 W] 


(5) 


where  PVap  is  the  solute's  vapor  pressure  at  the  given  temperature,  and  v  solid,  its  molecular 
volume  in  the  condensed  (solid)  phase.  The  numerator  in  (5)  is  the  ratio  of  actual  to  ideal- 
gas-predicted  solubility  (enhancement  factor).  The  denominator  is  called  the  Poynting 
correction  factor.  The  enhancement  factor  changes  by  several  orders  of  magnitude  as  the 


1 
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fluid  is  compressed  isothermally  from  a  low-density  gas  to  supercritical  densities.  In 
contrast,  the  Poynting  correction  is  quite  insensitive  to  pressure  (e.g.,  at  300K,  and  for  a 
typical  vsoiid  =  1.66  x  10'28  m3,  the  Poynting  factor  changes  from  1  to  1.49  in  going  from 
1  to  100  bar).  Thus,  the  reciprocal  of  the  solute's  fugacity  coefficient  is  directly 
proportional  to  the  enhancement  in  solubility  with  respect  to  ideal  gas  conditions  at  the 
given  temperature  and  pressure.  The  solute's  fugacity  coefficient  and  its  reciprocal  are  both 
plotted  versus  bulk  density  for  the  attractive  mixture  in  Figure  6.  The  extreme  sensitivity  to 
changes  in  density  is  obvious. 

To  examine  the  relationship  between  solubility  enhancement  and  microstructure,  we 
use  a  variable  radial  cutoff,  R,  in  the  calculation  of  the  integrals  that  lead  to  the  solute’s 
chemical  potential  (or,  equivalently,  to  its  fugacity).  Beyond  R,  the  mixture  is  assumed  to 
be  uniform  (gjj  =1).  Within  R,  pair  correlation  functions  from  integral  equation  calculations 
were  used.  This  calculation  provides  a  direct  measure  of  the  extent  of  the  region 
surrounding  the  solute  molecules  that  contributes  to  that  species'  fugacity  coefficient,  and 
hence  to  its  solubility.  The  radial  cutoff  required  in  the  evaluation  of  the  integrals  in  order  to 
obtain  a  quantity  equal  to  99.5,  99,  and  98%  of  the  true  (bulk)  fugacity  coefficient  is 
plotted  in  Figure  7.  There  is  a  noticeable  dependence  of  the  cutoff  upon  bulk  density.  As 
the  critical  density  is  approached,  larger  cutoffs  are  required.  It  is  obviously  the 
microstructure  (3<R*<5.5,  or  llA<R<20A)  that  determines  the  solute's  fugacity 
coefficient.  While  other  thermodynamic  quantities  in  dilute  supercritical  solutions  (such  as 
the  solute's  partial  molar  properties)  are  determined  by  long-ranged  correlations  in  the 
compressible  near-critical  region,  we  find  that  it  is  only  the  microstructure  around  the  solute 
(10-20  A)  that  affects  solvation,  and  hence  solubility. 

The  above-described  integral  equation  calculations  have  allowed  the  investigation  of 
the  spatial  distribution  of  solvent  and  solute  molecules  in  two  binary  Lennard-Jones 
mixtures:  one,  attractive,  with  potential  parameters  chosen  to  represent  naphthalene  in 
carbon  dioxide;  the  other,  repulsive,  with  potential  parameters  chosen  to  represent  neon  in 
xenon.  In  particular,  integral  equations  provide  a  direct  link  between  microstructure  (radial 
distribution  functions)  and  thermodynamics  (chemical  potentials  and  solubility). 
Accordingly,  we  have  placed  emphasis  in  understanding  how  differences  between  local  and 
bulk  conditions  around  solute  molecules  influence  bulk  thermodynamic  behavior.  As  with 
other  types  of  mixtures,  the  microstructure  is  important  in  determining  the  solute's 
solubility  in  a  given  solvent,  and  hence  we  also  investigated  its  relationship  to  fugacity 
coefficients.  What  is  unique  to  supercritical  mixtures  is  the  coexistence  of  widely  different 
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relevant  length  scales,  with  microstructure  controlling  solubility,  and  long-ranged, 
critically-driven  fluctuations  controlling  the  temperature  and  pressure  derivatives  of 
solubility.  In  the  attractive  mixture,  the  microstructure  is  found  to  be  solvent-rich  with 
respect  to  the  bulk.  This  enrichment  persists  even  when  the  solvent  density  is  averaged 
over  three  solvation  shells,  and  is  particularly  pronounced  in  the  density  range  0.5  <  p/pc  < 
0.8.  At  higher  densities,  packing  constraints  become  progressively  important,  the  pair 
correlation  function  gradually  acquires  liquid-like  oscillatory  character,  and  density 
enhancement  in  one  solvation  shell  (peaks)  are  canceled  by  depletions  in  the  successive  one 
(valleys).  While  the  magnitude  of  the  enhancement  in  the  first  solvation  shell  is  unrelated  to 
proximity  to  criticality,  the  persistence  length  of  this  density  perturbation  is  directly  affected 
by  the  growth  of  the  correlation  length.  Thus,  the  local  density,  when  averaged  over  two  or 
three  solvation  shells,  is  the  result  of  the  subtle  interplay  between  specific  (short-range)  and 
universal  (long-ranged)  effects.  The  solvation  region  around  the  attractive  solute  was  found 
to  be  a  sensitive  function  of  bulk  density,  and  varied  between  1 1  and  20  A  (3  and  5.5 
solvent  diameters),  with  the  maximum  occurring  near  the  solvent's  critical  density.  The 
repulsive  microstructure  was  found  to  be  solvent-lean.  Correlation  holes  appear  as  a  result 
of  competitive  attractions.  Deviations  between  local  solvent  density  (averaged  over  two 
solvation  shells)  and  bulk  conditions  are  particularly  pronounced  near  the  solvent's  critical 
density. 


Molecular  Dynamics  Simulations 

In  order  to  investigate  the  dynamic,  time-dependent  behavior  that  corresponds  to  the 
above-discussed  microstructural  features,  we  conducted  a  series  of  molecular  dynamics 
simulations  on  a  model  attractive  Lennard-Jones  mixture.  Potential  parameters  are  listed  in 
Table  III. 


Table  III.  Lennard-Jones  Parameters  for  Pyrene  in  CO2 


Pyrene 


in  CO2 

ij 

Oij  (A) 

Oij/022 

£ij/k  (K) 

£ij/£22 

solute-solute 

11 

7.140 

1.882 

662.8 

2.942 

solvent-solvent 

22 

3.794 

1.000 

225.3 

1.000 

solute-solvent 

12 

5.467 

1.440 

386.4 

1.715 

A  mass  ratio  of  mj/m2  =  4.597  was  used  to  simulate  the  pyrene-carbon  dioxide  system. 
The  simulations  were  carried  out  in  the  canonical  ensemble  (constant  volume,  V;  number  of 
molecules,  N;  and  temperature,  T).  A  time-step  of  0.004  02^m2  /  e2  ,  which  is  equivalent 
to  7.355  x  10' 15  sec.,  and  a  5th-order  predictor-corrector  scheme  with  an  automated  Verlet 
neighbor  list  (Chialvo  and  Debenedetti,  1991)  were  used  to  integrate  the  equations  of 
motion.  A  sample  size  and  solute  mole  fraction  of  4000  and  0.0025,  respectively,  were 
used  throughout  the  study  (10  solute,  3990  solvent  molecules).  The  simulation  length  at  the 
three  state  points  studied  so  far  was  200,000  steps  (1471  picoseconds).  These  are 
exceptionally  long  simulations,  a  requirement  imposed  by  the  need  to  gather  good  solute- 
solute  statistics  in  a  dilute  mixture.  Three  simulations  have  been  carried  out  so  far,  at 
densities  (  p*  =  po\  )  0.15,  0.31,  and  0.6,  and  at  a  temperature  (T*  =  kT/£2)  1.4.  This 
corresponds  to  a  solvent  reduced  density  range  of  0.48  to  1.93,  and  a  solvent  reduced 
temperature  of  1.07. 

In  the  course  of  each  simulation  we  studied  solute-solute  aggregation  statistics.  To 
this  end,  two  solute  molecules  were  considered  "bound"  at  any  instant  if  their  centers  were 
within  a  distance  not  exceeding  1.122  CT]  (location  of  the  solute-solute  potential  well).  The 
algorithm  of  Sevick  et  al.  (1988)  was  used  to  determine  connectivities,  and  hence  the 
instantaneous  distribution  of  solute-solute  aggregates.  Results  from  these  calculations  are 
summarized  in  Table  IV. 
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Table  IV.  Aggregate  Statistics  and  Collision  Dynamics  for  Pyrene-CC>2 


Observable 

p*  =0.15 

p*  =0.31 

p*  =  0.60 

10' 32  coll,  rate  (l1  s*1)^) 

1.16±0.27 

1. 8H0.45 

3.03H.24 

1011  k  (1  s'1  mol'1)^) 

14.813.5 

5.4H1.37 

2.4210.99 

%  monomers 

97.27+0.36 

98.1910.04 

98.0810.48 

%  dimers 

2.7710.59 

2.1010.28 

1.7910.46 

%  trimers 

0.16+0.10 

- 

0.17710.158 

dimer  lifetime  (ps) 

0.75810.045 

0.78U0.052 

0.79910.076 

trimer  lifetime  (ps) 

0.49410.292 

- 

0.42010.103 

(a)  Solute-solute  collision  rate,  r. 

(b)  Rate  constant  for  solute-solute  collisions,  k  =  r/pj2 

It  can  be  seen  that  the  solute  is  present  overwhelmingly  in  unassociated  form,  and 
that  this  behavior  is  unaffected  by  density.  Likewise,  the  average  lifetime  of  a  dimer  is  quite 
insensitive  to  density  variations.  Figure  8  shows  the  distribution  of  dimer  lifetimes  for  the 
three  densities  investigated  so  far.  It  can  be  seen  that  the  location  of  the  main  peak  is 
insensitive  to  density,  although  the  lifetime  distribution  changes  shape,  becoming  broader 
and  with  an  incipient  tendency  towards  bimodality  at  high  density. 

The  dynamics  of  the  solvation  shell  around  each  solute  molecule  were  studied  by 
computing  an  autocorrelation  function,  <D(t),  defined  as 


<aV(Q&V(0)> 

<8N(0)8N(0)> 


(6) 
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where 


T 

<  N(r)  >=  lim  f  N(r,t)dt 


(7) 


and  N(r,t)  is  the  instantaneous  number  of  solvent  molecules  located  within  a  sphere  of 
radius  r  centered  around  a  solute  molecule.  Thus,  8N  is  the  instantaneous  fluctuation  of 
N(r,t)  about  <N(r)>.  Figure  9  shows  the  behavior  of  <1>  at  high  density,  with  r  taken  to  be 
the  distance  to  the  first  minimum  of  the  solute-solvent  pair  correlation  function  (roughly 
2.55oi2).  The  autocorrelation  curve  is  averaged  over  the  ten  solute  molecules.  A 
convenient  measure  of  the  relaxation  time  for  the  solvation  shell  is  a  time  t,  defined  so  that 
d>(t)  =  0  (shortest  time  for  the  vanishing  of  the  autocorrelation).  The  density  dependence 
of  T  is  shown  in  Figure  10.  The  error  bars  show  the  range  of  values  corresponding  to 
individual  solute  molecules.  Two  features  of  T  are  noteworthy.  First,  the  characteristic  time 
for  the  decay  of  solvation  shell  density  fluctuations  is  quite  large  (of  the  order  tens  of 
picoseconds).  Secondly,  the  mean  1  shows  a  maximum  at  the  solvent's  critical  density. 
Although  variations  among  individual  solute  molecules  are  quite  large,  this  suggests  the 
type  of  dynamic  sluggishness  that  is  normally  associated  with  critical  slowing  down.  The 
origin  of  this  effect  is  the  growth  of  the  correlation  length  and  the  resulting  increased 
cooperativity.  Clearly,  longer  simulations  must  be  conducted  in  order  to  verify  this  trend. 

From  the  molecular  dynamics  study  we  conclude  that  the  very  high  solute-solute 
peaks  do  not  translate  to  solute-solute  aggregates;  that  the  solute  is  present  overwhelmingly 
in  unassociated  form;  that  this  behavior  is  insensitive  to  density;  that  the  typical  lifetime  of  a 
solute  dimer  (0.8  ps)  is  insensitive  to  a  four-fold  bulk  density  change;  and  that  the 
characteristic  time  for  density  fluctuations  within  the  solvation  shell  of  a  solute  molecule  is 
quite  long,  of  the  order  of  tens  of  picoseconds. 

4.  Current  and  Future  Work 

We  are  extending  the  molecular  dynamics  studies  to  higher  temperatures.  We  plan 
to  initiate  molecular  dynamics  simulations  of  molecular  fluids,  and  to  begin  the  planned 
collaboration  with  Dr.  James  Gord,  of  the  Systems  Research  Laboratory  in  Dayton,  who 
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will  be  investigating  supercritical  solvation  dynamics  experimentally  using  an  Ultrafast 
Laser  System.  Initially,  we  will  study  solvation  and  rotation  dynamics  in  the  benzene- 
carbon  dioxide  system  at  high  solute  dilution  and  supercritical  conditions  with  respect  to  the 
solvent  carbon  dioxide). 
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Figure  2:  Density  dependence  of  the  solute-solvent  pair  correlation  function  for  the 
naphthalene-carbon  dioxide  system  atT*  =  1.4  and  yi  =  10‘9. 
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Fieure  4:  Density  dependence  of  the  solute-solvent  pair  correlation  function  for  the  neon- 
xenon  system  at  T*  =  1.4  and  yi  =  lO9. 
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Figure  6:  Density  dependence  of  the  solute's  fugacity  coefficient  and  its  reciprocal 
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Figure  8:  Distribution  of  dimer  lifetimes  for  the  pyrene-carbon  dioxide  system  at  T*  =  1  4 
and  yi  =  0.0025. 
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N=4000,  N^IO,  T'=1.4,  p“= 0.60 


Figure  9:  Behavior  of  the  solvation  shell  autocorrelation  for  the  pyrene-carbon  dioxide 
system  at  T*  =  1.4,  p*  =  0.6,  and  yi  =  0.0025. 
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